Non-equilibrium dynamics of stochastic point processes with refractoriness 
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Stochastic point processes with refractoriness appear frequently in the quantitative analysis of 
physical and biological systems, such as the generation of action potentials by nerve cells, the 
release and reuptake of vesicles at a synapse, and the counting of particles by detector devices. 
Here we present an extension of renewal theory to describe ensembles of point processes with time 
varying input. This is made possible by a representation in terms of occupation numbers of two 
states: Active and refractory The dynamics of these occupation numbers follows a distributed delay 
differential equation. In particular, our theory enables us to uncover the effect of refractoriness on 
the time-dependent rate of an ensemble of encoding point processes in response to modulation of the 
input. We present exact solutions that demonstrate generic features, such as stochastic transients 
and oscillations in the step response as well as resonances, phase jumps and frequency doubling 
in the transfer of periodic signals. We show that a large class of renewal processes can indeed be 
regarded as special cases of the model we analyze. Hence our approach represents a widely applicable 
framework to define and analyze non-stationary renewal processes. 

PACS numbers: 02.50.Ey, 87.19.11, 29.40.-n 



I. INTRODUCTION 

Point processes are stochastic models for time series of 
discrete events: a particle passes through an apparatus, 
a photon hits a detector, or a neuron emits an action 
potential [1, 2]. As diverse as these examples are, they 
share three basic features that need to enter a statisti- 
cal description and which are illustrated in Fig. 1. The 
first feature is refractoriness. Technical devices to detect 
point events typically cannot discriminate events in ar- 
bitrarily short succession. This is addressed as the dead- 
time of the detector [3, 4] . The process of vesicle release 
and transmitter recycling in the synaptic cleft is of sim- 
ilar nature [5]. Upon the arrival of an action potential 
at the synapse, a vesicle might fuse with the membrane 
and release its contents into the synaptic cleft. Subse- 
quently the vesicle is reassembled for future signaling, 
but it is available only after a certain delay, equivalent 
to a refractory signalling component. In neurons, refrac- 
toriness can be the result of the interplay of many cellu- 
lar mechanisms, and possibly also of network effects [6]. 
In case of cortical neurons, which are driven to produce 
an action potential mainly by fluctuations of the input 
currents [7], refractoriness can model the time it takes 
to depolarize the membrane from a hyper-polarized level 
that follows the action potential into a range in which 
action potentials can be initiated by fluctuations. Gen- 
erally, refractoriness can be described as a duration d 
for which the component cannot be recruited to generate 
another event. In Fig. 1 it is illustrated as a delay line. 
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After the refractory time is elapsed, the component reen- 
ters the pool of active components that can generate an 
event. The existence of such a pool is the second com- 
mon property of the examples. Each component process 
of the ensemble can be either active or refractory. So 
an ensemble of neurons, vesicles, or detectors, can be 
treated in terms of the occupation of two states, "active" 
and "refractory", as depicted in Fig. 1, where A{t) <G [0, 1] 
describes the fraction of components which are active at 
time t and 1 — A{t) is the fraction of components that 
are currently refractory. The third feature is the stochas- 
tic nature of event generation. The time of arrival of a 
particle at a detector, the fluctuation of the membrane 
potential of a neuron that exceeds the threshold for ac- 
tion potential initiation, and the release of a vesicle into 
the synaptic cleft can under many conditions be assumed 
to happen stochastically. Given an independent transi- 
tion density of \{t) per time interval, event generation 
follows an inhomogeneous Poisson process, as indicated 
in Fig. 1. In the example of a detector, \(t) corresponds 
to the actual rate of incoming particles, and we will call 
it the input rate in the following. We distinguish two 
models of systems which share the properties described 
above: If the refractoriness has a fixed duration we obtain 
the well-known Poisson process with dead-time (PPD). If 
the duration is drawn randomly from a specified distribu- 
tion we call the model the Poisson process with random 
dead-time (PPRD). 

In the following we describe an extension of renewal 
theory for ensembles of point processes with time varying 
input. For stationary input rate, many previous publi- 
cations have investigated the statistics of the PPD [8- 
11]. In case of slowly varying input rates, expressions for 
mean and variance of detector counts have been derived 
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Figure 1. Scheme of the ensemble description of the Pois- 
son process with refractoriness: Active component processes 
produce events with rate X(t) and remain refractory for the 
duration d, illustrated by the delay line. After the dead-time 
they become active again. The fraction of active component 
processes is given by A(t). 



[12], and recently a method was proposed to correct for 
clustered input events [13]. PPDs with non-stationary 
input rate have been employed as a model for the signal 
transduction in the auditory nerve [14]. A sudden switch 
of the input rate was found to induce strong transients 
of the ensemble output rate. These reflect a transiently 
perturbed equilibrium of the occupation numbers in our 
two-state model, and we quantitatively analyze this case 
for the PPD. Response transients to rapidly changing 
input, in fact, explain the relation between neural refrac- 
toriness and neural precision [15]. Furthermore, periodic 
input profiles are known to be distorted by refractori- 
ness [16]. Here we derive the mapping of periodic input 
to output in the steady state and uncover the impact 
of refractoriness on the transmission. Interacting popu- 
lations of refractory neurons have been studied in [17]. 
This approach, in contrast to ours, neglects the effects 
of refractoriness on short time scales due to temporal 
coarse-graining of the population dynamics. 

From a more abstract perspective, the PPD is a very 
simple example of a point process that exhibits stochastic 
transients, which are not shared by the ordinary Pois- 
son process. Besides its many applications, the PPD 
therefore is a prototype system to study non-equilibrium 
phenomena in point process dynamics. Generally, non- 
stationary point processes can be defined by two different 
models: by rescaling time [18-20] or by time-dependent 
parameters of the hazard. The drawback of the former 
method is that the transformation from operational to 
real time distorts the inter-event intervals, such that, for 
example, a constant refractory period is not maintained. 
An example how a time-dependent hazard function can 
be derived from an underlying neuron model with time- 
dependent input is given in [21]. Our approach differs 
with regard to the choice of the hazard function, which 
enables rigorous analysis of the dynamics of the process. 

To analytically investigate non-equilibrium phenomena 
in ensembles of renewal process, a typical approach is to 
use a partial differential equation (PDE) for the proba- 
bility density of the ages of the components (time since 
the last event) [6]. In Section II we derive the two-state 
representation of the PPD from the dynamics of the age 
density. We present analytical solutions of the popula- 



tion dynamics for the response to a step change in the 
input rate in Section III, and to periodic input rate pro- 
files in Section IV. Finally, in Section V, we generalize 
our results to random refractoriness. We compute the 
effective hazard function of the resulting inhomogeneous 
renewal process, connecting it to the framework of re- 
newal theory. For the PPRD with gamma-distributed 
dead-times, as applied recently to model neural activity 
[22] , we show how the dynamics in terms of a distributed 
delay differential equation can be reduced to a system 
of ordinary (non-delay) differential equations. Again we 
study the transient response of an ensemble of processes 
to a step-like change in the input rate and the transmis- 
sion of periodic input. We observe that both distributed 
and fixed refractoriness lead to qualitatively similar dy- 
namical properties. At last we identify the class of re- 
newal processes that can be represented as a PPRD. As 
it turns out, this covers a wide range of renewal processes. 



II. DYNAMICS OF AN ENSEMBLE OF PPDS 

Point processes can be defined by a hazard function 
h(t, H t ) d = lim -P [event in [t, t + e]\U t ], (1) 

e— 5-0 e 

which is the conditional rate of the process to generate 
an event at time t, given the history of event times % t 
up until t. A process is a renewal process [1] if the haz- 
ard function depends only on the time r since the last 
event (age) instead of the whole history % t . This can be 
generalized to the inhomogeneous renewal process which, 
additionally, allows for an explicit time dependence of the 
hazard function h{t,W t ) = h(t,r). 

Here we consider an ensemble of point processes de- 
fined by the hazard function 

h(t,r) = \(t)6(r - d), (2) 

where 9{t) = {1 for t > 0, else} denotes the Heaviside 
function, d > is called dead-time, A(i) > is the time- 
dependent input rate and r > is the age of the compo- 
nent process. This is an inhomogeneous renewal process, 
which is known as the Poisson process with dead-time 
(PPD). The state of an ensemble of such processes can 
be described by the time-dependent probability density 
of ages a(t, r), for which a partial differential equation is 
known [6] 

d d 

— a(t,r) = -—a(t,r) - h(t,r)a(t,T). (3) 

Solutions must conserve probability, which manifests it- 
self in the boundary condition a(t,0) — v(t), with the 
event rate of the ensemble 

u(t) = f h(t, r)a(t, t) dr = \{t)A(t). (4) 
Jo 
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In the second step we inserted (2) and introduced the 
active fraction of component processes with age r > d 



A(t) 



def / 
Jd 



a(t, t) dr. 



(5) 



For t < d, Eq. (3) simplifies to d t a(t,r) = — <9 T a(i, t), 
implying 



a(t + u, u) = a(t, 0) = v(t) Vue[0,d). 



(6) 



Since a(t, r) is normalized we obtain with the boundary 
condition and (6), 

1= / a(t,T)dT= / v(s)ds + A(t). (7) 

JO Jt~d 

This equation is the starting point of the analysis of in- 
teracting populations of refractory neurons in [17]. Dif- 
ferentiation of (7) by t yields 



d 
dt 



A(t) = X(t - d)A(t - d) - X(t)A(t) 



(8) 



which is a linear delay differential equation (DDE) with 
time-dependent coefficients. Its forward solution for in- 
put A(t) is uniquely defined given A(t) on an interval of 
length d [23]. However, not all solutions of (8) can be 
interpreted physically, since by differentiation of (7) ad- 
ditive constants are lost. Only if the initial trajectory 
satisfies (7) , Eq. (8) determines the time evolution of the 
ensemble. With (4), the time-dependent output rate v(t) 
follows. Note that only in the case of the "pure" Pois- 
son process with d = we obtain v(t) = X(t), because 
A{t) = 1 by (7). 

Eq. (8) represents a more accessible description of the 
process in terms of the occupation of the active and the 
refractory state (see Fig. 1 and Section I) compared to 
the dynamics of the probability density of ages (3) . This 
description is feasible because of the particular nature of 
the hazard function of the PPD (2) . In the following we 
will consider specific solutions of (8). 



III. SOLUTIONS FOR A STEP INPUT 

If A(t) = A is constant, given the occupation A(t) = 
u(t) on the first interval t e [— d, 0] with u : \—d 1 0] — > 
[0, 1], solutions to (8) are known in integral form [23] 



A{t) = u(0)g(i) + / Xu(s - d) g(t - s) ds (9) 
Jo 

for t > 0, where we introduced the fundamental solution 
g(t). It obeys g(t) — {0 for t < 0, 1 for t = 0} and solves 
(8) for t > 0. As we will show, here g is in fact the shifted 
and scaled auto-correlation function R of the process. 

The inter-event interval density of the stationary PPD 
is f(t) = \0{t - d)e- x(t - d \ For t > d, the integral equa- 
tion 



A(t) = (/*A)(t), 



(10) 



is equivalent to the delay differential equation (8), which 
can be proven by differentiation with respect to t (* de- 
notes the convolution). The auto-correlation function [1] 



(ii) 



k=0 



with f* k (t) = (/*(*-!) * f)(t) for k > 1 and f*° d = S(t), 
solves (10) for t > d, and hence is a solution of (8) in 
that domain. We find for k > 1 that 



f* k (t) = \ k {t-kd)- 



sfc-l. 



-A(t- 



kd ^e(t-kd)/(k-i)\. (12) 



Given the initial trajectory g(t) for t < 0, solving (8) by 
variation of constants for t <G [0, d] yields g(t) — X~ 1 f(t + 
d) = \~ 1 R(t + d). Then due to uniqueness of the solution 
it holds for all t > that 



g(t) = X^R(t + d) 



(13) 



We apply these results to compute the response of A(t) 
if the input rate is switched from A to A at t = 0, given 
the process was in equilibrium for t < 0. Eq. (7) deter- 
mines this equilibrium to A(t) = ao = (l + Ao<i) _1 , t < 0. 
In this case, the step change in X(t) enters (9) as 

A stev {t) = u(0)g(t) + X(s- d) u(s - d) g(t - s) ds, 

Ao 

(14) 



for t > 0. We insert (13) to obtain 



Istcp 



(t) 



— R{t + d) + — 



[ R(t 
Jo 



+ d — s) ds 



ao 
~X 
a-oXo 
~X~ 



R(t + d) 



AqOo 



1 



1 



R(t + d) 



X \ X 
(I + {Xq 1 - X- 1 ) R(t + d)) 



(15) 



where we used (7), which holds for g(t) = X~ 1 R(t + d). 
Fig. 2 shows this analytical solution compared to direct 
numerical simulation of an ensemble of PPDs upon a step 
change of the input rate X(t) at t = 0. The output rate 
displays a marked transient, which increases with the 
dead-time d and exhibits oscillations of frequency 1/d. 



IV. TRANSMISSION OF PERIODIC INPUT 

We now investigate an ensemble of PPDs with an input 
rate X(t) £ R that is periodic. If T is its period, we obtain 
the Fourier series X(t) = Efel-oo A fc e lfcut , with uj = ?f 
and At e C. Then the steady state solution for the active 
fraction A(t) of the PPD is also periodic in T, so it can 



be expressed as A(t) — J2T=- 
Inserted into (7) we obtain 



a k e lku)t with a k e 



1= £ ^ k qi +k e l{l+k)u,t + J2 ake 



ikut 



(16) 
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Figure 2. Transients upon step change of the input rate X(t) at t = 0. Exact analytical result (15) (solid lines) and simulation 
of the ensemble rate of 10 10 processes (crosses). Parameters: d[s] : 0.02, 0.05, 0.08 (light gray, mid gray, dark gray) A: 
A = ((5Hz)- 1 - d)-\ A = ((10Hz)- 1 - d)"\ B: A = ((lOHz)" 1 - c() _1 , A = ((5Hz)- 1 - d)' 1 . 



where for k ^ 



1 _ p -iktud 
Jkojt (]£ = „ikuit def 



t-d 



ikuj 



-e = q k e 



ikuji 



and qo = d. Since the Fourier basis functions {e lkujt , k G 
Z} are mutually orthogonal, we can separate (16) for 
different k. This yields the infinite dimensional linear 
system of equations 



h,o = Qk A ' a *=-' + a k, k G Z. 



(17) 



The ensemble averaged output rate of the PPD defined 
in (4) then follows 
with the spectrum 



in (4) then follows as v(t) = X(t)A(t) = £fct-oo (3 k e ikult 



Pk = J! a k-l^l = % 1 (<^fe,0 - a k) ■ 



(18) 



£=- 



where we used (17). This given, we replace the a k by (3 k 
in (17) to obtain 



/3fe = A fe - A k-mPi, kez. 



(19) 



i=- 



This relation shows how different frequencies of the out- 
put rate are coupled by a convolution with the input 
spectrum. Note that inverting (19) yields the spectrum 
of the time-dependent input rate A(t) given the spectrum 
of the output rate signal v(t). 

Let us now consider the special case of a cosine- 
modulated input 

X(t) = Xq + ecos(ix>t), 

which we obtain with = {0 for |fc| > 1, | for k G 
{1,-1}, A for k = 0}, A > e > 0. Then for k € N, 
(17) becomes a so-called three-term recurrence relation 
[24] of the form = a n+ i + x n a n + y n a n -i with x n = 



(g r 7 1 + Ao)(2/e) and y n = 1. This relation has two linearly 
independent solutions. The unique minimal solution is 
convergent and can be obtained from the continued frac- 
tion r„_i = —y n / {x n +i"n) in a robust manner [24] using 
the relation r n — a n+ i/a n , n > 0: Setting rjy = for 
some N € N one computes (r n )o< n <N backwards and 
increases N until tq does not change within the required 
tolerance. Inserting a\ — tqchq into (17) for k = we 
solve for «o to obtain ao = (1 + c?(Ao + e3ff(ro))) _1 (here 
3? denotes the real part). The remaining a k follow recur- 
sively from a k +i = a k r k and a- k = ct\, since Ait) G R. 
The spectrum of the output rate is then given by (18). 
Fig. 3A shows the output rate v(t) for different input rate 
modulation frequencies / = lo/(2tt). Fig. 3C,D display 
the amplitude and phase of the three lowest harmonics 
of the output rate v(t) as a function of /. The time 
averaged emission rate (/3q) depends on the modulation 
frequency. It is maximized slightly below the character- 
istic frequencies / = k/d. This is due to the oscillation of 
A{t), which is almost in phase at these frequencies and 
hence cooperates with the oscillatory hazard rate A(i) 
to enhance the emission (see Fig. 3D). Interestingly, the 
first and second (^2) harmonic of v(t) display max- 
ima at different /. At a particular modulation frequency 
/ ~ l/(2d) the amplitude of the second harmonic (^2) is 
larger than the first harmonic (/?i), so that the ensem- 
ble activity is effectively modulated with twice the input 
frequency (see Fig. 3A (a) and Fig. 3C): the ensemble 
performs a frequency doubling. Fig. 3B shows the maxi- 
mum over one period of the output rate trajectory. These 
maxima are dominated by the maxima of the amplitude 
of the first harmonic. In particular, low frequency in- 
put signals are transmitted to the output with strong 
distortion and reduced intensity, because the fraction of 
non-refractory processes, A(t), is in anti-phase (Fig. 3D) 
to X(t) and hence suppresses the output rate's modula- 
tion. This is in contrast to the common view that the 
PPD transmits slow signals more reliably than the Pois- 
son process [6[. Note that only if the driving frequency 
/ = n/d, n G N is an integer multiple of the inverse 
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dead-time then A(t) = (1 + Aqc?) -1 is constant in time 
and the output rate is proportional to A(t) without any 
distortion. 

V. RANDOM DEAD-TIME 

For detector devices as well as for neurons, a fixed 
dead-time might be a somewhat restricted model. Here 
we consider the PPRD as described in the introduction. 
Upon generation of each event, the PPRD draws an in- 
dependent and identically distributed random dead-time 
with the probability density function (PDF) p for the 
duration of which it remains silent . The PPRD is still a 
renewal process, since it has no further dependencies on 
the event history beyond the time since the last event. 
As in the case of a fixed dead-time in Section II, the fol- 
lowing analysis of the PPRD is based on the conservation 
of the total number of processes in an ensemble. Inactive 
components must have generated an event at some time 
in the past, which leads to the normalization condition 

/t />oo 
A(t')X(t') / p{x)dxdt' . (20) 
-oo Jt-t' 

This equation can be seen as the generalization of the 
normalization condition (7) , from which the DDE (8) fol- 
lows by differentiation, to the case of random dead-time. 
Analogously the distributed DDE 

—A(t) = -X(t)A(t) + / p(x)\(t-x)A(t-x)dx (21) 

follows from (20) by differentiation with respect to t. 
Eq. (21) describes the time-evolution of the occupation of 
the active state for an ensemble of general PPRDs. Obvi- 
ously, the dynamics of the PPD (8) is recovered from (21) 
in case of the localized density p(x) — 5(x—d). In the rest 
of this section, we will derive the hazard function h(t, r) 
of the PPRD, consider the case of gamma-distributed 
dead time and the associated step response, generalize 
the transmission of periodic input to random dead-time, 
and finally identify the class of renewal processes that 
can be represented by the PPRD. 

For a given density p(x) of the dead-time it is not ob- 
vious what the hazard function of the PPRD is. In order 
to relate the PPRD to renewal theory we compute its 
time-dependent hazard function (1) here. Let 

Q(t, t) d = E [6(t - a;) | last event at t - t] 

denote the probability of the process to be active at time 
t, given the last event occured at t — r, where x is the 
random dead-time and E denotes the expectation value 
with respect to x. The hazard function is then h(t, r) = 
\(t)Q(t,T). With 

Q(t, t) = P [x < t I last event at t — r] 
= 1 — P [x > t | ev. at t — t n no ev. in (t — r, t)] 
= 1 - P [x > t] /P [no ev. in (t -r,t)\ ev. at t - r] 



we obtain 

h(t, r) = A(t) (1 - J-(r)/E [F(t, r\x)}) , (22) 

where 

/oo 
p(x)dx 

is the survivor function of the dead-time distribution and 
F(t, t\x) = cx P (-(9(t - x) [ \{t')dt') 

Jt-T+X 

is the survivor function of a PPD with dead-time x. In 
case of constant X{t) = A we further have 

E [F(t, t\x)} = e- XT f e Xx p(x)dx + T(t). 
Jo 

The hazard function (22) is shown for constant X(t) in 
Fig. 4A for the special case described below. Eq. (22) was 
applied to generate realizations of the PPRD for Fig. 4B. 

For gamma-distributed dead-time (8) can be trans- 
formed into a system of ordinary differential equations. 
We exemplify the application of (21) for gamma dis- 
tributed dead-times with parameters neN and (3 € K + , 

p(x) = K n (x), (23) 
Kn(ar) = /3 n+1 x n e- 0x /n\, (24) 

with E[x] = (n + l)//3. The time course of the rate can 
be obtained from (21). Introducing 

/t 
Kk(t — x)v(x) dx 
-oo 

for < k < n and b n+1 (t) d = A(t) and exploiting the 
relation 

^K k (x) = 6(k - l)^/c fc _ 1 (a;) - f)K k (x) 
dx 

for < k < n enables to replace the integral in (21) by a 
closed system of ordinary differential equations 

(-b n+1 (t)X(t)+b n (t) k = n+l 
jb k (t) = I |86 fc _i(t) -Ph(t) 1 < k < n (25) 

[/3b n+1 (t)\(t) - (3b (t) k = 0. 

For constant X(t) = X this can be written as -^b(t) = 

M x b(t), b(t) e M™+ 2 . Hence given the initial state 6(0) 
the solution unfolds to 

b{t) = cxp(M A t) • &(0) . (26) 

With v(t) = v = (A" 1 +EIX})- 1 the equilibrium state fol- 
lows: Setting the temporal derivatives to in (25) yields 
bk = v for < k < n, and b n+1 = A = l—i/E[x]. The rate 
response to a switch from A to A at t = is thus given 
by (26) where 6(0) is the equilibrium state for A . A nu- 
merical simulation of the process with gamma distributed 
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Figure 3. Transmission of cosine-modulated input. A, B: Theoretical result (solid lines) and simulation of an ensemble of 10 10 
processes (crosses). A: Steady-state rate v(t) for different modulation frequencies /, with fd: 0.42, 0.85, 1.0, 1.4 (a,b,c,d). 
B: !Anax = max(z/(i)) for different /. Here d[s): 0.02, 0.05, 0.08 (light gray, mid gray, dark gray) C,D: Amplitude (C) and 
phase (D) of harmonics k £ {0, . . . , 3} of A(t) (19) (top) and v(t) (18) (bottom) as a function of modulation frequency /. 
Grayscale denotes order of harmonics k : 0, 1, 2, 3 (black, dark gray, mid gray, light gray), d — 80ms. Other parameters in A-D: 
\{t) = A (l + 0.9cos(2tt/£)), Ao = (i^~ 1 - d)"\ v = 10Hz. 



refractoriness (23) with hazard function (22) and the cor- 
responding analytical solution (26) upon a step change of 
X(t) are shown in Fig. 4. The simulation of the process 
was done via rejection [25] and averaged over indepen- 
dent runs. The spread of dead-times (Fig. 4A) does not 
qualitatively change the shape of the response transient 
(Fig. 4B). 

Analogous to Section IV, we consider the case of peri- 
odic input. We insert the Fourier series of X(t) and A(t) 
into (20) and obtain the same relation of their spectra 
(16) as for a single dead-time with the altered coefficients 



<lk 



roc pc 

= / e~ ikuy 
JO Jv 



p{x)dxdy. 



(27) 



As is easily seen, by inserting the localized dead-time 
PDF p(x) = 5(x — d) the original are recovered. Hence 
all results of Section IV also hold for the PPRD, but with 
the coefficients (27). In particular we would like to em- 
phasize the validity of the general input-output mapping 
(19) for arbitrarily distributed dead-time. 

Let us now investigate which class of renewal pro- 
cesses can be represented by the PPRD. We start with 
an arbitrary renewal process with inter-event interval 
/£R+, defined by its PDF l(x). Let E > be an in- 
dependent, exponentially distributed interval with PDF 
e(x) = Ae~ Xx , and let R be the random dead-time with 
PDF p(x). For / to be a realization of a PPRD it must 



hold for some p and A > that 

I = R + E^L = p-k€^L = f>€ 

=> p = A _1 (s + A)i => p = A _1 £ _1 [st] + l 
1 / d 



p{x) 



A V dx 



i(x) + i(0) )+i(x), (28) 



where * decorates a function which was transformed by 
the Laplace transform C, and s denotes the Laplace vari- 
able. The renewal process defined by i can be represented 
by a PPRD if p is a PDF. Let us call the hazard function 
of the renewal process h(x), and the survivor function 
F(x) = cxp(— J Q h(x')dx'), which obey t(x) = h(x)F(x) 
[1]. Assume that l(x) is differentiable. Since expression 
(28) is always normalized, in order for it to define a suit- 
able PDF we only have to require p(x) > for all x. 
possibly in the sense of distributions. This translates 
into 



\- 1 (h'(x)-h 2 (x))+h(x)>0. 
In case h(x) > 0, this can be written as 
, , x h'(x) , 



(29) 



(30) 



If, in addition, the hazard and its derivative are bounded 
in the sense that h(x) < oo and h'(x) > — oo, there exists 
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a A > such that (30) is fulfilled. These conditions are 
indeed met by a large class of renewal processes. 

For example, the gamma-process which has random 
inter-event intervals with PDF t(x) — tt r (x) (23) with 
parameters r, /? € R, r > 1, j3 > is by (28) equivalent 
to the PPRD with p{x) = n r -i{x) and A = /3, but other 
choices of A are also possible. This illustrates the well 
known fact that the inter-event intervals of gamma pro- 
cesses with integer shape parameter n can be considered 
as the concatenation of n exponentially distributed in- 
tervals. In neuroscience the gamma-process is frequently 
used to model stationary time series of action potential 
emissions of nerve cells. To describe adaptation phenom- 
ena, a time-dependence of the parameters of the haz- 
ard function was introduced in [21]. Identification of the 
gamma-process with a PPRD entails the alternative to 
generalize the gamma process to time-dependent rates by 
varying the input rate of the PPRD. Similarly, the log- 
normal process can be represented as a PPRD. We define 
its inter-event interval as x — £A, where £ is a unit-less 
random number and A gives the time-scale. Let £ be 
distributed according to the log-normal PDF 

^7^ cxp l 

for £ > 0, ?7(0) = 0, with unit-less parameters fi, a. Then 
x is distributed according to l(x) — A~ 1 t]( x /a). Accord- 
ing to (28) the process can be represented by any of the 
PPRDs with 

^ = 4.4-^(1 + ^)) 

A > A^V- 2 exp (-1 - \i + a 2 ) , 

where the lower bound on A is due to the requirement 
p{x) > 0. For these and other renewal processes for which 
a PPRD representation exists, non-equilibrium dynamics 
can be studied on the basis of (21). 



VI. DISCUSSION 

In this paper, we consider the effect of refractoriness 
on the output of an encoding point process in case of 
arbitrary time-dependent input signals. Such point pro- 
cesses, for example, are used to model the generation 
of action potentials by nerve cells, the release and reup- 
take of vesicles into the synaptic cleft, or the detection of 
particles by technical devices. We describe ensembles of 
these stochastic processes by the occupation numbers of 
two states: active and refractory. The active components 
behave as inhomogeneous Poisson processes, but after an 
event is produced the component is silent for the duration 
of the dead-time, it is caught in a delay line. We derive a 
distributed delay differential equation that describes the 
dynamics in the general case of a randomly distributed 
dead-time. 



Due to the simpler dynamics in case of a fixed dead- 
time, we first elaborate properties of the PPD. For sta- 
tionary input rate, we solve the dynamics of the en- 
semble in a way that sheds light on the connection be- 
tween the fundamental solution of the DDE and the auto- 
correlation function of the point process. This relation 
is employed to express the time-dependent ensemble rate 
(output) for a step-change of the hazard rate (input). 
The resulting output rate displays stochastic transients 
and oscillations with a periodicity given by the dead- 
time. Such transients might enable nerve cells to respond 
reliably to rapid changes in the input currents [15, 26]. 
For periodically modulated input rate, we demonstrate 
how the spectrum of the steady-state periodic output rate 
results from the linear coupling between harmonics. In 
the particular case of cosine-modulated input signals only 
adjacent harmonics are coupled. This nearest-neighbor 
interaction is rigorously solved using the theory of three- 
term-recurrence relations and continued fractions [24]. 

Our analytic result explains frequency doubling, the 
emergence of higher harmonics and the dependence of 
the time averaged population activity on the modulation 
frequency. In particular, slow frequency components of 
the input are attenuated and distorted in the popula- 
tion rate, which is in contrast to the claim that the PPD 
transmits slow frequency signals more reliably than the 
Poisson process [6]. 

In case of periodic input modulation, the output spec- 
trum contains all harmonics of the fundamental fre- 
quency of the input. This might be related to a psy- 
chophysical phenomenon called "missing fundamental il- 
lusion" [27, 28]: Being presented an auditory stimulus 
which consists of several harmonics of a fundamental fre- 
quency, but in which the fundamental frequency itself is 
missing, subjects nonetheless perceive the fundamental 
frequency as if it was contained in the stimulus spec- 
trum. By considering neurons in the auditory system as 
PPDs whose hazard rate is modulated by the auditory 
stimulus, our theory explains how the lowest harmonic 
is recovered in the population activity of the neurons. 
Conversely, our results can be applied to infer input rate 
profiles from the count rate of detectors with dead-time, 
in particular in the case of periodic input, for which (19) 
applies. 

For the more general case of a random, arbitrarily dis- 
tributed dead-time, we show how the DDE generalizes 
to a distributed DDE. By suitable choice of the distri- 
bution of the dead-time, non-equilibrium dynamics of a 
large class of renewal processes can be described. For in- 
teger gamma-distributed dead-time we demonstrate how 
the distributed DDE transforms into a coupled system of 
finitely many ordinary differential equations, which could 
also be implemented as a multi-state Markov system [22] . 
Regarding the output rate transient upon a step change 
of the input and the transmission of periodic inputs, we 
find that the qualitative behavior of the system is very 
similar to the PPD. In conclusion, we present a canon- 
ical model for non-stationary renewal processes, as well 
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Figure 4. PPD with random dead-time with mean 80ms, shape parameter n: 10, 50 (black, gray). A: density of dead- 
times p(r)/ max(p(r)) (23) (dotted lines) and hazard function h(r)/Xo (22) for X(t) = Xq (solid lines). B: Transients upon 
step change of the input rate X(t) at t = 0. Theoretical result from (26) (solid lines) and simulation of an ensemble of 10 6 
processes with hazard function (22) (crosses) averaged over 225 trials. The error bars denote the standard deviation over trials. 
A = ((5Hz)- 1 - d)-\ X = ((lOHz)- 1 - d)- 1 . 



as the analytical methods to describe ensembles thereof. funded by BMBF Grant 01GQ0420 to BCCN Freiburg 
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